In [60]:
import numpy as np
import astropy
from pearce.mocks import cat_dict
from pearce.mocks.assembias_models.table_utils import compute_prim_haloprop_bins
import h5py
from os import path
In [61]:
from matplotlib import pyplot as plt
%matplotlib inline
import seaborn as sns
sns.set()
from itertools import cycle
In [62]:
config_fname = '/home/users/swmclau2/Git/pearce/bin/mcmc/nh_gg_sham_hsab_mcmc_config.yaml'
In [63]:
import yaml
with open(config_fname, 'r') as ymlfile:
cfg = yaml.load(ymlfile)['data']
In [64]:
sim_cfg = cfg['sim']
obs_cfg = cfg['obs']
cov_cfg = cfg['cov']
In [65]:
#sim_cfg['simname'] = 'ds_14_b_sub'
#sim_cfg['halo_property'] = 'halo_vmax@mpeak'
sim_cfg['gal_property_fname'] = '/scratch/users/swmclau2/smf_dr72bright34_m7_lowm.dat'
In [66]:
cat = cat_dict[sim_cfg['simname']](**sim_cfg['sim_hps']) # construct the specified catalog!
# TODO logspace
r_bins = obs_cfg['rbins']
obs = obs_cfg['obs']
if type(obs) is str:
obs = [obs]
meas_cov_fname = cov_cfg['meas_cov_fname']
emu_cov_fname = cov_cfg['emu_cov_fname']
if type(emu_cov_fname) is str:
emu_cov_fname = [emu_cov_fname]
assert len(obs) == len(emu_cov_fname), "Emu cov not same length as obs!"
n_bins = len(r_bins)-1
data = np.zeros((len(obs)*n_bins))
assert path.isfile(meas_cov_fname), "Invalid meas cov file specified"
try:
cov = np.loadtxt(meas_cov_fname)
except ValueError:
cov = np.load(meas_cov_fname)
assert cov.shape == (len(obs)*n_bins, len(obs)*n_bins), "Invalid meas cov shape."
In [80]:
cat.filenames[-1]
Out[80]:
In [82]:
from halotools.sim_manager import RockstarHlistReader
In [86]:
reader = RockstarHlistReader(cat.filenames[-1], cat.columns_to_keep, cat.cache_filenames[-1],\
cat.simname, 'rockstar', 0.0, cat.version_name, cat.Lbox, cat.pmass,\
overwrite=False, header_char = '#')
In [87]:
reader.read_halocat(cat.columns_to_convert)
In [89]:
reader.halo_table
Out[89]:
In [67]:
#cat.load(sim_cfg['scale_factor'], **sim_cfg['sim_hps'])
#cat.populate()# will generate a mock for us to overwrite
gal_property = np.loadtxt(sim_cfg['gal_property_fname'])
halo_property_name = sim_cfg['halo_property']
min_ptcl = sim_cfg.get('min_ptcl', 200)
nd = float(sim_cfg['nd'])
scatter = float(sim_cfg['scatter'])
In [91]:
cat.h
Out[91]:
In [90]:
from AbundanceMatching import *
af = AbundanceFunction(gal_property[:,0], gal_property[:,1], sim_cfg['af_hyps'], faint_end_first = sim_cfg['reverse'])
remainder = af.deconvolute(scatter, 20)
# apply min mass
halo_table = reader.halo_table#cat.halocat.halo_table#[cat.halocat.halo_table['halo_mvir']>min_ptcl*cat.pmass]
nd_halos = calc_number_densities(halo_table[halo_property_name], cat.Lbox) #don't think this matters which one i choose here
catalog_w_nan = af.match(nd_halos, scatter)
n_obj_needed = int(nd*(cat.Lbox**3))
catalog = halo_table[~np.isnan(catalog_w_nan)]
sort_idxs = np.argsort(catalog)
In [69]:
np.all(cat.halocat.halo_table['halo_upid']==-1)
Out[69]:
In [92]:
np.sum(halo_table['halo_upid'] == -1)
Out[92]:
In [93]:
final_catalog = halo_table[~np.isnan(catalog_w_nan)][sort_idxs[:n_obj_needed]]
#final_catalog['x'] = final_catalog['halo_x']
#final_catalog['y'] = final_catalog['halo_y']
#final_catalog['z'] = final_catalog['halo_z']
#final_catalog['halo_upid'] = -1
# FYI cursed.
#cat.model.mock.galaxy_table = final_catalog
# TODO save sham hod "truth"
In [94]:
final_catalog.colnames
Out[94]:
In [95]:
from halotools.mock_observables import *
In [96]:
x, y, z = [final_catalog[c] for c in ['halo_x', 'halo_y', 'halo_z']]
pos = return_xyz_formatted_array(x, y, z, period=cat.Lbox)
In [97]:
xi_all = tpcf(pos / cat.h, r_bins, period=cat.Lbox / cat.h, num_threads=4,
estimator='Landy-Szalay')
In [98]:
xi_all
Out[98]:
In [100]:
plt.plot(xi_all)
plt.yscale('log')
In [102]:
import h5py
f = h5py.File('/home/users/swmclau2/scratch/PearceMCMC/pearce_mcmc_nh_gg_sham_hsab.hdf5', 'r')
In [103]:
cov = f['cov'].value
In [105]:
np.savetxt('/home/users/swmclau2/Git/pearce/bin/shams/xigg_cov_mcmc.npy', cov)
In [ ]: